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Abstract 

Background: Tumor necrosis factor (TNF) is a widely studied cytol<ine (ligand) tliat induces proinflammatory 
signaling and regulates myriad cellular processes. In major illnesses, such as rheumatoid arthritis and certain 
cancers, the expression of TNF is elevated. Despite much progress in the field, the targeted regulation of TNF 
response for therapeutic benefits remains suboptimal. Here, to effectively regulate the proinflammatory response 
induced by TNF, a systems biology approach was adopted. 

Results: We developed a computational model to investigate the temporal activations of MAP kinase (p38), nuclear 
factor (NF)-kB, and the kinetics of 3 groups of genes, defined by early, intermediate and late phases, in murine 
embryonic fibroblast (MEF) and 3T3 cells. To identify a crucial target that suppresses, and not abolishes, 
proinflammatory genes, the model was tested in several in silico knock out (KO) conditions. Among the candidate 
molecules tested, in silico RIPl KO effectively regulated all groups of proinflammatory genes (early, middle and late). 
To validate this result, we experimentally inhibited TNF signaling in MEF and 3T3 cells with RIPl inhibitor, Necrostatin-1 
(Nec-1), and investigated 10 genes {116, Nfl<bia, Jun, TnfaipS, Ccl7, Vcaml, CxcllO, Mmp3, Mmpl3, Enpp2) belonging to the 
3 major groups of upregulated genes. As predicted by the model, all measured genes were significantly impaired. 

Conclusions: Our results demonstrate that Nec-1 modulates TNF-induced proinflammatory response, and may 
potentially be used as a therapeutic target for inflammatory diseases such as rheumatoid arthritis and osteoarthritis. 
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Introduction 

The tumor necrosis factor (TNF), first termed in 1962 [1], 
was initially known for its ability to induce programmed 
cell death or apoptosis. As a result, throughout the years, 
the TNF has been intensely investigated for its anticancer 
property [2]. Today, this cytokine is central to the regula- 
tion of myriad important cellular processes such as prolif- 
eration, differentiation, growth, and the immune response. 

TNF binds to two types of outer membrane bound re- 
ceptors on target cells, TNFRl and TNFR2, and triggers 
the cell survival and proinflammatory NF-kB and MAP 
kinases activations [3]. In addition, the TNFRl induces 
intracellular cell death pathways via caspases after intern- 
alization through endocytosis. It is, therefore, conceivable 
that the dysregulation of the TNF signaling process will 



* Correspondence: kumar@ttck.keio.ac.jp 

^Institute for Advanced Biosciences, Keio University, 14-1 Baba-cho, Tsuruoka, 
Japan 

^Systems Biology Program, Graduate School of Media and Governance, Keio 
University, Fujisawa, Endo 5322, Japan 



misbalance proinflammatory and/or apoptotic responses. 
Notably, the chronic aberration in the baseline levels of 
TNF in human circulatory system has been attributed 
to the pathogenesis of numerous diseases, including 
rheumatoid arthritis, osteoporosis, sepsis and cancer [4,5]. 

The vast majority of TNF related biological processes 
are initiated by the death-domain (DD) containing 
TNFRl, which is also called TNFRSFIA. Unlike TNFR2, 
TNFRl is present in almost all cell types in humans. Upon 
TNF binding, TNFRl trimerizes, and its intracellular DD 
recruits TRADD, which then creates a platform for RIPl 
and TRAF2 to collectively form the receptor-signaling 
complex I. Cellular inhibitor of apoptosis proteins (cIAP)- 
1 and -2 bind to complex I and, consequently, together 
with K63-linked ubiquitin chains, modify RIPl and 
TRAF2 [6]. This creates docking sites for an E3 ligase or 
linear ubiquitin chain assembly complex (LUBAC) consist- 
ing of heme-oxidized IRP2 ubiquitin ligase- 1 (HOIL-1), 
HOIL-1 -interacting protein (HOIP), and SHANK- 
associated RH domain interacting protein (SHARPIN). 



o 
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Subsequently, the activation of TAKl and the ubiquitina- 
tion of NEMO (or IKKy), a subunit of IKK complex, lead 
to cell survival or proinflammatory response through 
NF-kB and MAP kinases activations. Other TRAP super- 
family members (TRAPS and 6) are also known to play a 
role in the NP-kB and MAP kinases activations [7,8]. 

On the other hand, for the apoptotic pathways, clathrin, 
AP-2 and Dyn first mediate receptor internalization. 
Receptor-signaling complex I becomes modified, and dis- 
sociates from TNPRl, allowing PADD and caspase-8 to 
form complex II. Within complex II, caspase-8 becomes 
activated to induce extrinsic apoptosis through caspase-3 
activation. Alternatively, caspase-8 activates caspase-7, and 
eventually, the cleavage of Bid to tBid in the mitochondria 
activates caspase-9 via cathepsin D. This induces the in- 
trinsic apoptosis through caspase-3 activation. 

Due to its ability to signal numerous cellular processes 
via the survival and death pathways, the TNPRl signaling 
research has received immense attention over the years, 
especially on understanding the downstream signaling cas- 
cades to regulate and control proinflammatory diseases 
and cancer. Despite numerous studies, the control of pro- 
inflammatory diseases through therapeutic treatments, 
where TNP is over-expressed, remains suboptimal. Por ex- 
ample, biologic response modifiers or biologies, such as 
Etanercept and Infliximab, are TNP decoy receptors or 
antibodies that suppress TNPRl signaling through compe- 
tition for TNP. Although these drugs have shown success- 
ful downregulation of inflammation in many cases, they 
can immuno-compromise patients to secondary infections 
such as tuberculosis [9], or have been ineffective in a sub- 
stantial number of administered patients [10]. 

To find alternatives, there have been major efforts on se- 
lectively suppressing the intracellular signaling of TNPRl. 
Por example, genetic knockouts (KOs) of TRAPs and 
TRADD acting on the proinflammatory pathways have 
been investigated [7,8,11]. However, the experimental out- 
comes, so far, have not been optimistic. In TRAP2 KO, 
there is compensatory activation of NP-kB through TRAPS 
[7] or TRAP6 [8], and vice-versa. On the other hand, 
TRADD KO almost completely abolishes NP-kB activation 
[11], which is not desirable for the general survivability of 
cells. Thus, a systemic approach where the propagation of 
signal transduction to all known branching pathways dur- 
ing target intervention should be monitored. This will 
allow the elucidation of effective target candidate(s) that 
overcomes and balances the deficiencies of current 
investigations. 

In this paper, we adopted a systems biology approach to 
study TNPRl signaling dynamics. Pirstly, we developed a 
computational model of TNP-induced proinflammatory 
response leading to NP-kB, MAP kinase activations, and 
three groups of gene expressions (classified according to 
their temporal profiles [12]). The model is based on the 



perturbation-response approach [13-16], which has been 
successfully used to elucidate novel signaling features and 
behaviors in Toll-like receptor-4 [17,18], -3 [19], and TNP- 
related apoptosis -inducing ligand (TRAIL) signaling [20]. 
Secondly, the TNPRl model parameters were selected to fit 
the temporal activation profiles of NP-kB and MAP kinase 
p38 for fibroblast cell type in several available conditions 
(wildtype [7], TRAP2 KO [7], TRAPS KO [7], TRAP2/ 
TRAPS double KO (DKO) [7], TRAP6 KO [8], TRADD 
KO [11] and RIPl KO [21]). Using the resultant TNPRl 
model with robust parameters, we performed simulations 
of multiple in silico KOs to determine an optimal target 
that suppresses, but not abolishes, proinflammatory genes. 
Pinally, to validate the modeling results, we performed ex- 
periments measuring various key proinflammatory gene ex- 
pressions in MEP and 3T3 cells for TNP stimulation. 
Overall, our study presents evidence that systems biology 
research can be useful to elucidate important target(s) to 
suppress proinflammatory diseases such as rheumatoid 
arthritis and osteoarthritis. 

Results 

TNFR1 signaling topology and model 

To develop a computational model of proinflammatory 
TNPRl signaling dynamics, we first require the known 
signal transduction pathways. We curated the KEGG data- 
base, and performed literature survey of the latest TNP re- 
search. After carefully considering several sources, we 
were able to propose a signaling topology mainly by com- 
bining the knowledge from KEGG, Palschlehner et al. 
(2012) and Wertz et al. (2010) [6,22] (Pigure 1). 

Next, to simulate TNP-induced dynamics of NP-kB 
and MAPK activations using the topology, we devel- 
oped a dynamic model based on perturbation-response 
approach (Materials and Methods), using COPASI 
simulation platform [23]. Unlike common biochemical 
reaction models [24,2S], the perturbation-response ap- 
proach does not require detailed knowledge of all signal- 
ing species and their reaction kinetics. This is because it 
analyses the response waves of signal transduction instead 
of individual reaction kinetics [13- IS, 17-20]. The response 
waves can be approximated using linear response rules 
{Response Rules ^ Additional file 1: Pigure SI) combined 
with the law of mass conservation, and this approach has 
been previously used to successfully model the TLRs and 
TRAIL signaling pathways [17-20]. 

Briefly, each reaction in the model is represented by 
a first-order response equation with activation or de- 
activation term. The activation term generally refers to 
protein binding, transformation, complex formation, 
phosphorylation and transcription. The deactivation 
term refers to protein unbinding, dephosphorylation 
and negative regulation such as mRNA decay through 
microRNA regulation. 
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(See figure on previous page.) 

Figure 1 Schematic of TNFR1 signaling of cell survival/proinflammatory and apoptosis pathways. Upon TNF receptor activation, 
complexes I (survival pathway) and II (apoptosis) are formed. Complex I subsequently activates transcription factors, such as activator protein 
(AP)-l and NF-kB through MAP kinases and IKK complex, respectively, which subsequently bind to promoter regions of genes to induce 
numerous proinflammatory genes. 

V J 



Simulating TNF-induced NF-kB and MAP kinase dynamics 

The parameters of the initial model (rate constants, or the 
elements of Jacobian matrix /, Materials and Methods) 
were estimated by fitting the simulation profiles with ex- 
perimental profiles of signaling molecules where data is 
available. We obtained published semi-quantitative experi- 
mental profiles of iKBa phosphorylation (NF-kB activa- 
tion) and p38 (MAP kinase) activation in wildtype and 
various genetically mutant MEFs generally treated with 
10 ng/mL of TNF (Figure 2 A, Additional file 1: Figure S2 
and Table SI) [7,8,11,21]. (Note that the kinetics of other 
MAP kinases, JNK and ERK, were also similar to p38 
[7,8,11,21]. Thus, we used p38 as a representative MAP 
kinase for our investigation). 

The parameter values were selected by using Genetic 
Algorithm [26] module in COPASI software [23] to fit the 
experimental profiles (Figure 2A, WT). Following, we per- 
formed sensitivity analysis (Materials and Methods) of the 



model parameters and found them to be robust to a 
small degree of uncertainty to their values (Additional 
file 1: Table S2). As a further validity of the parameter 
values, we tested the wildtype model in other condi- 
tions, namely TRAF2 KO, TRAF5 KO, TRAF2/5 double 
KO, TRAF6 KO, RIPl KO and TRADD KO (Figure 2B). 
(Note that in silico KOs were generated from the wild- 
type model by setting the activation parameter value of 
the KO molecule to null). 

Remarkably, we were able to obtain a single set of model 
parameters (Table 1, reactions 1-29 and see Additional file 
2 for the TNFRl model A in SBML format), which could 
be used to simulate the semi- quantitative profiles of iKBa 
phosphorylation and p38 kinase activation in multiple ex- 
perimental conditions. In wildtype, TRAF2 KO, TRAF5 
KO and TRAF6 KO, the iKBa phosphorylation and p38 
kinase activation reach peak values around 15 min and 
gradually decay at 30 min. Notably, TRAF6 KO shows 
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Figure 2 Experimental and simulated profiles of IkBq and p38 activations in wildtype and mutant conditions. (A) Experimental profiles, 
MEFs were generally treated with 10 ng/mL of TNF, and (B) simulated profiles of IkBq (top panels) and p38 (bottom panels) activations up to 
30 min in wildtype {W\), TRAF2 KO, TRAF5 KO, TRAF2mAF5 double KO aRAF2/5 DKO), TRAF6 KO, TRADD KO and up to 15 min in RIPl KO. Note: 
p38 experimental profiles are available only for MMJ, TRAF6 KO and TRADD KO. Experimental details and data are found in references [7,8,1 1,21]. 
ImageJ was used to estimate the intensities of the activation dynamics (Additional file 1: Table SI) for each molecule in each condition relative to 
wildtype peak activation values found in Additional file 1: Figure S2. IkBq phosphorylation refers to NF-kB activation throughout the text. 
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(1-12): in-sHico knock-out conditions are performed by setting parameter values (/c,) to 0 for targeted reactions in TRADD KO (1), clAP1/2 KO (2), TRAF2 KO (3), TRAF5 KO (4) TRAF2/5 DKO (3,4), RIP KO (5), TRAF6 KO 
(6), TAK1 complex KO (7), LUBAC KO (8), SHARPIN KO (9), IkBq KO (10), MKK3/6 KO (11) and p38 KO (12). (13) Kinetics of Group II mRNA transcription and decay processes were refitted after adding feedback (without 
feedback: /C35 = 7e-3, /C40 = 1.2e-5). Bold italic fonts (reactions 48-65) indicate additional feedback activation pathways required for group III continuous activation. * indicates the multiplication sign. 
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enhanced iKBa phosphorylation and p38 kinase activation 
due to Signaling Flux Redistribution {Response Rule 5, 
Additional file 1: Figure SI) [18]. In the remaining condi- 
tions, the activation levels of both molecules are very weak 
(RIPl KO and TRAF2/5 DKO) or absent (TRADD KO). 

It is noteworthy that although there have been previous 
models on TNF signaling [24,27,28], to our knowledge, 
this is the first time a single model of TNF signaling with 
fixed parameter values recapitulates the proinflammatory 
signaling dynamics in multiple experimental conditions. 

To compare our linear response model (TNFRl model 
A) simulations with other models that contain more de- 
tafled descriptions of IKK [28] and MAPK [29] signaling, 
using higher order terms and Michaelis-Menten type 
kinetics, we developed an alternative TNFRl model B in- 
corporating the relevant reaction details (Additional file 1: 
Table S3). Notably, the simulations of TNFRl models A 
and B show very similar dynamics for a fixed amount of 
TNF perturbation (Additional file 1: Figure S3). Thus, we 
concur that our linear response model can be appropri- 
ately used for further investigations. 

Simulating distinct TNF-induced gene expression patterns 

Next, we extended the TNFRl model (we will now simply 
call TNFRl model A as TNFRl model) to simulate down- 
stream proinflammatory gene expression dynamics. Re- 
cently, time-series high throughput microarray and 
quantitative real time PGR experiments on TNF simulated 
mouse 3T3 fibroblasts cells have revealed 3 distinct groups 
of upregulated gene expression patterns, with possibly cor- 
responding distinct biological roles [12,30]. The groups 
were labeled into "early I", "intermediate or middle 11" and 
"late III" response, according to their time to reach peak 
expressions between 0.5-1, 2-3, and 6-12 h, respectively, 
after TNF stimulation (Figure 3A) [12,30,31]. Here, we ex- 
tended the TNFRl model to simulate the temporal profiles 
of the 3 groups of gene expressions. 

According to our modeling approach, the time to peak 
activation can be controUed by reaction parameter 
values and/or the number of signaling intermediates 
[15,17-20]. Briefly, decreasing (increasing) the activation 
or transcription parameter value wiU show lower 
(sharper) gradients of formation part of the expression 
profiles. Alternatively, decreasing (increasing) the deacti- 
vation or decay parameter value will show lower 
(sharper) gradients of depletion part of expression pro- 
files {Response Rule i. Additional file 1: Figure SI). In 
addition, inserting intermediary reactions between tran- 
scription process and gene induction will increase delay 
for gene expression dynamics {Response rule 2, Add- 
itional file 1: Figure SI). The intermediates can represent 
the complexities of transcription process involving the 
pre-initiation, initiation, promoter clearance, elongation 
and termination [32], or post- transcriptional processes 



such as messenger RNA editing and splicing. Using this 
approach, the TNFRl model was extended to simulate 
the temporal dynamics of groups I, II and III genes. 
Note that the response rules (Additional file 1: Figure SI) 
are used to modify an initial signaling topology only 
after aU parameter space has been exhaustively searched, 
and a reasonable model fit is unable to be achieved [20]. 

Previous investigations on the 3 groups of genes have 
indicated distinct mechanisms for the differential dy- 
namical response [12,30]. Hao and Baltimore have found 
lesser presence of AU Rich Element (ARE) region on the 
3'UTR of group III genes, targeted by microRNAs and 
ARE-binding proteins (such as tristetraprolin) that en- 
hance RNA decay processes. Hence, it was postulated as 
one possible reason for the lower decay response of 
group III genes compared with genes from groups I and 

II [12]. More recently, by studying the kinetics of pre- 
mRNA and mRNA, Hao and Baltimore observed delays in 
splicing of groups II and III genes compared to group I 
genes. The differential delays were suggested as another 
biological mechanism for the distinct gene profiles [30] . 

In our extended model, we, therefore, considered both 
mechanisms to reproduce the temporal profiles of the 3 
groups of genes. Notably, our simulations of pre-mRNA 
and mRNA for aU groups of genes matched the data of 
Hao and Baltimore for the first 60 min (Additional file 1: 
Figure S4). However, subsequently for 12 h, although the 
simulations of groups I and II genes were recapitulated, 
group III simulation was poor (Figure 3B, blue line). 
Specifically, reducing the parameter value for the decay 
term representing lower miRNA and ARE-binding pro- 
teins regulating decay processes {Response rule i. Add- 
itional file 1: Figure SI), and adding intermediates 
{Response rule 2, Additional file 1: Figure SI) to provide 
delays in RNA splicing in our model were not sufficient 
to produce the continuous activation of group III genes 
(Figure 3C, cyan dotted line). 

To overcome the shortfall in the model simulations, 
we hypothesized that novel activation or transcription 
term(s) (positive feedback) may be present to provide 
additional flux for the continuous increase in group 

III expressions {Response rule 4, Additional file 1: 
Figure SI). This could result from secondary post- 
transcriptional/translational mechanisms through i) 
autocrine signaling such as IL-1 [33], IL-6 [34] or 
TGF-|3 [35] signaling (Figure 3D), or ii) cytosolic feed- 
back mechanisms specifically for group III genes [36] 
(Figure 3E). Thus, a novel feedback mechanism pre- 
dominantly affecting the transcription of group III genes 
was added to the TNFRl model (Table 1, equations 
30-65 and Additional file 2). 

The modified TNFRl model with feedback mecha- 
nisms to group III genes produced simulations that 
matched all 3 groups of gene expression profiles 
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Figure 3 Three distinct groups of TNF-activated genes. (A) Average expression profiles of genes in groups I (red, peal< at 0.5 in), II (green, 2h) 
and III (blue, 12h) in 3T3 fibroblasts stimulated with recombinant mouse TNF. Figure was reproduced from [12]. (B) Simulation profiles of the 3 
groups of genes using the initial TNFRl model. (C) Simulation of the modified TNFRl model with transcriptional delay and novel feedback 
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intermediates included in the updated TNFRl model (refer to Table 1, reactions 48-65), and Y refers to a novel transcription factor, such as 
interferon regulatory factor (IRF). 



(Figure 3A and C, solid lines). To scrutinize the feed- 
back mechanism, we re-monitored the simulation pro- 
file of NF-kB for 6 hours (Additional file 1: Figure S5). 
The resultant profile mimics the damped oscillatory 
dynamics of NF-kB previously observed in murine fi- 
broblasts [36]. Overall, these data suggest that low 
miRNA regulation and additional delay in RNA spli- 
cing are not sufficient to produce the continuous 
activation of group III genes, and that a novel tran- 
scription process, possibly through secondary post- 
transcriptional/translational autocrine signaling, such 
as IL-1 signaling or other novel feedback mechanisms 
that activate NF-kB, and not MAPK (Additional file 1: 
Figure S5), are required. 



Predicting key target for regulating proinflammatory 
response 

Now that the TNFRl model is able to successfully simu- 
late the three groups of upregulated genes in wildtype, 
we investigated the significance and effect of removing 
or suppressing key intracellular signaling molecules for 
controlling proinflammatory response, in silico. 

It is well known that TNFRl signaling is enhanced in 
proinflammatory diseases and cancer [1-4]. To investi- 
gate which known molecules would be potential target 
to regulate the cell survival or proinflammatory activity, 
we performed in silico KOs of all possible signaling mol- 
ecules within the TNFRl model. In total, we simulated 
groups I, II and III dynamic gene expressions in 12 
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(TRADD, cIAPl/2, TRAF2, TRAF5, TRAF6, RIPl, 
SHARPIN, LUBAC, TAKl complex (TAKl/TABl/2), 
iKBa, MKK3/6 and p38) KO conditions and compared 
with wildtype profiles (Additional file 1: Figure S6). 

Among the candidates, the removal of TAKl complex 
or RIPl produced the most noticeable downregulation 
of all 3 gene groups, which chiefly consist of well-known 
proinflammatory mediators (Figure 4). However, in 
TAKl complex KO, our simulations show almost no in- 
duction for group 1 genes. The substantial impairment 
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Figure 4 The effects of in sHico KOs on the expression profiles 
of the 3 groups of genes. Simulated expression profiles of group 1 
(A), group II (B), and group III (C) genes in 4 experimental 
conditions: wildtype {WT), IkBq KO, RIPl KO and TAKl complex KO 
for 12 hours using the modified TNFRl model (with feedback). 
Predictions of RIPl KO (orange curves) indicate suppression, but not 
abolishment, of all groups of gene expressions compared to 
wildtype (black curve). 



in gene expressions (> 90%) is usually detrimental to the 
general survivability of living cells, and this has been par- 
ticularly demonstrated in TAKl -deficient mice [37,38]. 
RIPl, on the other hand, showed about 50-70% impair- 
ment compared to wildtype peak expressions. Our simula- 
tions, therefore, suggest that RIPl is possibly a crucial 
single molecule target for controlling enhanced proinflam- 
matory response due to TNFRl signaling in proinflamma- 
tory disease conditions, such as in rheumatoid arthritis, 
without compromising the normal functioning of other 
cellular activities. 



Experimental inhibition of RIP1 downregulates 
proinflammatory genes in TNF stimulation 

To verify the predictions of TNFRl model simulations, 
we prepared corresponding MEF and BALB/3T3 ceUs 
treated with TNF in wfldtype and in RIPl suppression. 
Necrostatin-1 (Nec-1) was originally identified as a po- 
tent small molecule inhibitor of necroptosis or non- 
apoptotic cell death [39]. Further interests in Nec-1 led 
to its specificity towards the inhibition of RIPl [40]. Al- 
though Nec-1 has recently been extensively studied, its 
effect on the expressions of groups I, II and III genes in 
TNF stimulation remains largely unknown. Therefore, 
here, we used Nec-1 to suppress RIPl in vivo. 

To check the effect of ceU death by Nec-1, we compared 
MEF and BALB/3T3 cells treated with different doses of 
Nec-1 in the presence or absence of TNF (Additional file 
1: Figure S7). The data revealed that Nec-1 has no sub- 
stantial effect on cell death after 24 h incubation, and 
hence, could be tested for its efficacy on the 3 groups of 
genes. We next performed quantitative RT-PCR for a total 
of 10 genes: 7/6, TnfaipS, Jun, Nfl<bia (group I), Cc/7, 
Vcaml, CxcllO (group II), and Mmp3, MmplS, Enpp2 
(group III). We intentionally included key proinflamma- 
tory mediators, genes of matrix metalloproteinase {Mmp3, 
Mmpl3), which are known to degrade collagen in cartilage 
and thereby enhance rheumatoid arthritis and osteoarth- 
ritis progression [41-44]. 

A previous study has shown that 30 [iM of Nec-1 effect- 
ively inhibited RIPl kinase activity [41]. Therefore, we in- 
vestigated gene expressions for cells stimulated with 
10 ng/mL TNF, in the presence or absence of 30 [iM Nec- 
1 for a period of 10 hours with measurements made at 
least every hour (Figure 5). Remarkably, as predicted by 
the TNFRl model, RIPl inhibition by Nec-1 resulted in 
the suppression of aU 3 groups of genes. The effect of sup- 
pressing RIPl is significant for groups I and II genes in 
both MEF and BALB/3T3 cells, especially during the first 
2-3 hours after stimulation. For group III genes, Nec-1 
had more pronounced effect in MEF compared with 
BALB/3T3 cells. Overall, these results are consistent with 
the TNFRl model predictions that suppressing RIPl in 
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TNF stimulation significantly impairs the activation of all 
3 groups of genes. 

Discussion 

TNF is a crucial cytokine that regulates myriad vital cellular 
processes. However, its levels are enhanced in major proin- 
flammatory diseases. Here, to understand the TNF-induced 
proinflammatory signaling process, and to carefully regulate 



its dynamic response, a systems biology approach was 
adopted. We first developed a dynamic computational 
model using well-established publicly available experimen- 
tal data of NF-kB, MAP kinase p38, and the average profiles 
of 3 groups of 180 upregulated genes in mouse fibroblast 
cells. 

Despite the simplicity of using first-order response equa- 
tions to simulate the profiles of the intracellular molecules. 
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the computational model of TNFRl recapitulated the ex- 
perimental response in wildtype and several mutant condi- 
tions for NF-kB and p38 activations. This result is 
surprising, as we know that the innate immune response 
of TNF is highly complex. It is important to note here that 
there have been previous other computational efforts on 
NF-kB and MAPK signaling that had utilized detailed bio- 
chemical reactions modeling, to elucidate local properties 
of signal transduction, such as the ability of common mol- 
ecules to produce distinct feedback mechanisms to differ- 
ent stimuli [24,45,46]. In our work, however, we have 
shown that even a simpler representation of the signal 
transduction pathways, through first order response equa- 
tions and the law of mass conservation can reproduce ex- 
perimental dynamics. This strongly indicates the presence 
of simple organizing rules governing the deterministic 
population average signaling response [47-52] . 

Next, through the analyses of downstream temporal 
gene expression profiles, the model suggests the presence 
of additional novel post-transcriptional/translational pro- 
cesses that is required for the continuous activation of 
group III genes. This result is additional to previous postu- 
lations, which had indicated that the continuous activation 
is due to lesser ARE region for group III genes leading to a 
very low decay process [12], and due to the presence of 
differential delays in the RNA splicing process [30]. Our 
model suggests that, on top of these effects, a novel time- 
delayed secondary transcriptional mechanism is required. 

Literature survey indicates that the novel positive feed- 
back processes could be a result of autocrine signaling, ex- 
ample through IL-1 or IL-6, or derive from a still unknown 
intracellular feedback mechanisms regulating mainly the 
promoter regions of group III genes. For example, the role 
of interferon regulatory factor (IRF) family in inducing CclS 
or RANTES expression, which belongs to one of the group 
3 genes, is reported in a previous study [53], however, was 
not considered in the initial TNFRl model. It is, therefore, 
necessary to perform further experimental work to confirm 
and elucidate the exact mechanisms for the continuous ac- 
tivations of group III genes. 

On the other hand, for down-regulating TNF signal- 
ing, which is enhanced in several proinflammatory dis- 
eases and cancer, we performed the simulations for 12 in 
silico KOs of signaling molecules. The resultant simula- 
tions indicated that RIPl is a major regulator of the 3 
groups of upregulated gene expressions. To verify the re- 
sult, we performed experiments on MEF and BALB/3T3 
cells using Nec-1 as an inhibitor of RIPl. The measure- 
ment of 10 genes belonging to groups I {116, TnfaipS, 
Jun, Nfkbia), II (Cc/7, Vcaml, CxcllO) and III [MmpS, 
MmplSy Enpp2) all showed significant impairment with 
Nec-1 compared to wildtype. 

Most importantly, the expressions of key proinflam- 
matory genes such as 116, Vcaml, Ccl7, Mmp3, Mmpl3, 



enhanced in rheumatoid arthritis and osteoarthritis 
[41,44], were reduced. In particular are the levels of 
matrix metalloproteinase genes Mmp3, MmplS, which 
are known to directly affect type II collagen in bone car- 
tilages and degrade the extracellular matrix. Although 
recent therapeutics have been focusing on the specific 
regulations of MMPs [42-44,54], it remains to be seen 
what effect such treatments will have on other proin- 
flammatory or vital genes. 

In summary, our approach provides a systemic analysis 
of TNFRl signaling, and suggests Nec-1 is potentially an 
important therapeutic target for effectively regulating 
major proinflammatory mediators in chronic diseases 
where TNF is overexpressed. 

Materials and methods 

Computational model 

The model is based on perturbation-response approach 
[15,17-20]. The basic principle behind the approach is to 
induce a controlled perturbation of input reaction spe- 
cies of a system (TNFRl), and monitor the response of 
the activation/concentration levels of other output spe- 
cies (e.g. TAKl, p38, NF-kB, 7/6, etc.) from steady-state. 
To briefly explain the principle, let a stable network con- 
sisting of n species be perturbed from the reference 
steady-state. In general, the resultant changes in the con- 
centration of species are governed by the kinetic evolu- 
tion equation [13,14]: 



dt 



(1) 



where the corresponding vector form of equation 1 is 
^ = F{X) . F is a vector of any non-linear function in- 
cluding diffusion and reaction of the species vector X - 
(Xj^ X2, X^), which represents activated concentration 
levels of reaction species. The response to perturbation 
can be written by X = Xq + 5X, where Xq is the reference 
steady-state vector and 5X is the relative response from 
steady-states (^X^^q - 0). 

The generally non-linear kinetic evolution equation 1 
can be approximated or linearized by using Taylor series: 



MX 



+ . 



(2) 



As the general volume of perturbing substance is usu- 
ally very small (order of 1%) compared to the total vol- 
ume of cells that are perturbed [55], now consider a 
small perturbation around the steady-state in equation 2, 
in which higher-order terms become negligible and re- 
sult in the approximation of the first-order term. In vec- 
tor form ^ = ^ix^ \x=Xq ^-^ (note the change from 
partial derivative to total derivative of time), where the 
zeroth order term F{Xo) = 0 at the steady-state Xq 
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and the Jacobian matrix, or linear stability matrix, is 
J ~ ^ax"^ \x=Xq- elements of /, based on the initial ac- 
tivation topology, are chosen by fitting bX with corre- 
sponding experimental profiles. Hence, the amount of 
response (flux propagated) along a biological pathway 
can be approximated using first order mass- action re- 
sponse, i.e. ^ = J8X. That is, the basic principle so far 
suggests that the response rate of species in a mass- 
conserved system at an initial steady-state can be 
approximated by first order mass-action response equa- 
tions, given a small perturbation to one or more species. 

Note that Jacobian matrix elements (or response coef- 
ficients) can include not only reaction information, but 
also spatial information such as diffusion and transport 
mechanisms. Thus, each species in the perturbation- 
response model can represent a molecule, a different 
modified state of a molecule (e.g. ubiquitinated state) or 
a molecular process such as diffusion, endocytosis, etc. 
That is, each species in the biological network does not 
necessarily represent a specific molecular species. For il- 
lustration, in a pathway ^ X2 ^ X3 ^ ^ X5, Xi to 
Xs can each be a different species or the same species at 
different stages in signaling, for example, Xi being inter- 
nalized (becoming X2), transported to a different organ- 
elle (X3), ubiquitinated (X4) and become part of a 
protein complex (X5). 

The complete SBML version of TNFRl Models A & B 
are available in Additional file 2. 



Sensitivity analysis 

We performed a sensitivity analysis to test the robust- 
ness of the optimized model parameters using the 
COPASI sensitivities module with default values. The 
variation in the response of signaling molecules/steps, Xi 
(t), was analyzed when a small variation of each model 
parameter kj was applied. The response sensitivity coeffi- 
cient [56] of the molecule with regard to the par- 
ameter is defined by 



The obtained values, Rij are then scaled, to reflect 
the relative changes in response, such as a change of 
p% in the value of parameter kp results in a Rpj'P% 
change in the value of the peak activation of the 
molecule. The response sensitivity coefficients of p38, 
iKBa, and groups I, II and III genes were obtained at 
peak time {t = 15 min for p38 and IkB, 30 min, 2 h and 
12 h for groups I, II and III respectively, see Additional 
file 1: Table S2). 



Experiments 

Reagents and cell culture 

Recombinant mouse TNF was purchased from R&D sys- 
tems. Necrostatin-1 was purchased from Merck Millipore. 
3T3 cells were obtained from JCRB cell bank. 3T3 and 
MEF were grown in DMEM (Nissui Seiyaku Co.) contain- 
ing 10% calf serum, 100 U/mL of penicillin at 37°C in a 
5% CO2 humidified atmosphere. 

Evaluation of cell survival by 3-(4,5-dimethylthiazol-2-yl)- 
2,5-diphenyltetrazolium bromide (MTT) assay 

The sensitivity of cells to hyperosmotic stress was mea- 
sured with the MTT colorimetric assay in 96-well plates. 
Cells (2 X 10^) were inoculated in each well and incubated 
for 24 h. Thereafter, 50 \iL of MTT (2 mg/mL in PBS) was 
added to each well and the plates were incubated for a fur- 
ther 2 h. The resultant formazan was dissolved with 
100 \iL of dimethyl sulfoxide after aspiration of culture 
medium. Plates were placed on a plate shaker for 1 min 
and then read immediately at 570 nm using TEC AN mi- 
croplate reader with Magellan software (Mannedorf, 
Switzerland). 

Real-time PCR analysis 

Total cellular RNA was extracted from cells using the 
TRIzol reagent according to the manufacturer s instruc- 
tions (Invitrogen). One microgram of RNA was reverse- 
transcribed using a first-strand cDNA synthesis kit 
(ReverTra Acea; Toyobo). Quantitative real-time PCR 
was performed using SYBR premix Ex Taq (Takara) on 
the Applied Biosystems StepOnePlus™ according to 
the technical brochure of the company. RT-PCR 
primers designed in this study are listed in Additional 
file 1: Table S4. Quantitative measurements were deter- 
mined using the AACt method and expression of 
GAPDH was used as the internal control. Melt curve 
analyses of all real-time PCR products were performed 
and shown to produce the sole DNA duplex. 

Additional files 



Additional file 1: Figure SI. Response rules. Figure S2. Experimental 

raw data used for model fitting. Figure S3. Experimental vs. simulated 
profiles of IkBq and p38 activations in wildtype and mutant conditions 
using TNFRl model B. Figure S4. Simulation of pre-mRNA and mRNA ex- 
pression profiles of groups I, II and III genes. Figure S5. Simulation of NF- 
kB activation profiles with and without feedback mechanisms. Figure S6. 
The effects of in silico KOs on the expression profiles of groups I, II and III 
genes. Table SI. Estimation of the relative intensities of IkBq and p38 ac- 
tivation dynamics. Table S2. Sensitivity analysis of TNFRl model A. Table 
S3. TNFRl model B details. Table S4. List of primer sequences for RT-PCR. 

Additional file 2: TNFR1 models A & B in SBML format. 
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